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We predict a variety of composite quiescent and spinning two- and three-dimensional (2D and 3D) 
self-trapped modes in media with a repulsive nonlinearity whose local strength grows from center 
to periphery. These are 2D dipoles and quadrupoles, and 3D octupoles, as well as vortex-antivortex 
pairs and quadruplets. Unlike other multidimensional models, where such complex bound states 
either do not exist or are subject to strong instabilities, these modes are remarkably robust in the 
present setting. The results are obtained by means of numerical methods and analytically, using 
the Thomas-Fermi approximation. The predicted states may be realized in optical and matter-wave 
media with controllable cubic nonlinearities. 
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I. INTRODUCTION 

Spatial and spatiotemporal solitons in two- and three-dimensional (2D and 3D) settings are one of central research 
topics in photonics n-iii and physics of quantum gases H [5], as well as in other areas, such as nuclear matter [6] 
and the field theory [?]• This topic also provides a strong incentive for studies of soliton solutions and their stability 
in applied mathematics [8]. Typically, solitons result from the balance between diffraction and/or group-velocity 
dispersion, which drive linear defocusing of the fields, and their self-focusing induced by the attractive nonlinearity. 
While in one dimension (ID) this mechanism readily gives rise to stable solitons [Hll], a fundamental problem in the 
2D and 3D geometry is that the conventional Kerr (cubic) nonlinearity leads to the critical (in 2D) and supercritical 
(in 3D) collapse [9], hence the respective multidimensional solitons are unstable. 

Multipoles are a species of localized multidimensional modes which are of significant interest to nonlinear optics m 
and studies of matter waves in Bose-Einstein condensates (BEC) [11], but are also vulnerable to the collapse-induced 
instability. They appear as multiple-peak solitons with alternating signs of adjacent peaks, the simplest example being 
a dipole which consists of two lobes with opposite signs of the field in them. While, as mentioned above, multipoles 
are unstable in uniform Kerr media, possibilities were proposed to stabilize them, including the introduction of optical 
lattices [12], and the use of non-Kerr nmn] and nonlocal [15] nonlinearities. 

Related to these are schemes elaborated for the stabilization of multidimensional solitons in a more general context. 
One approach relies upon the use of more complex photonic nonlinearities, such as saturable [16] (which are available 
in photorefractive crystals HZ] and atomic vapors m), or cubic-quintic (CQ) [19] and quadratic-cubic [2] [20] com¬ 
binations of competing focusing and defocusing terms. Recently, the stability of 2D spatial solitons in optical media 
featuring the CQ [25] and quintic-septimal [26] nonlinear responses has been demonstrated experimentally. 

Nonlocal photonic nonlinearities also readily provide stabilization of multidimensional solitons against small per¬ 
turbations [21]. The solitons in nonlocal media are as well stable to thermal fluctuations [22] and partial loss of 
incoherence [23]. Peculiarities of the modulational instability driven by the nonlocal nonlinearity, which is the pre¬ 
cursor of the formation of solitons, were studied too [24] . 

The stabilization of 2D solitons may also be secured by “management” techniques, i.e., periodic alternation of the 
sign of the nonlinear term [27]. Still another possibility for the stabilization is provided by trapping potentials - in 
particular, periodic lattices, as has been predicted theoretically for various settings 13111 EH] and was demonstrated 
experimentally in photorefractive crystals [29] . 

A novel collapse-free scheme for the creation of fundamental and topologically structured D-dimensional solitons 
was proposed in Ref. [30] and further developed in a number of subsequent publications m - [36] . It is based on using 
the self-defocusing cubic nonlinearity with the local strength growing at r ^ oo at any rate exceeding r^, where r is 
the distance from the center. In reality, the local strength does not need to attain extremely large values, as the solitons 
supported by this setting are strongly localized, hence properties of the medium at distances essentially exceeding 
the transverse size of the soliton are unimportant m- In optics, the 2D version of the scheme can be implemented 
for spatial solitons by means of inhomogeneously doping the bulk waveguide by nonlinearity-enhancing resonant 
impurities m In fact, a uniform dopant density may be used, with an inhomogeneous distribution of detuning from 
the two-photon resonance imposed by an external field [30]. In 2D and 3D atomic BECs, the spatial modulation 
of the scattering length of two-body collisions, which determines the local nonlinearity strength in the mean-field 
approximation, can be imposed by means of the optically controlled spatially nonuniform Eeshbach resonance, as 
experimentally demonstrated in Ref. [38]. Alternatively, the required spatial profile of the nonlinearity can be 
“painted”, as an averaged one, by a fast moving laser beam [39]. Another experimentally demonstrated method 
makes use of the magnetically-controlled Eeshbach resonance in an appropriately shaped magnetic lattice, into which 
the condensate is loaded m- 

Eurthermore, the scheme for the creation of solitons in media with the spatially growing nonlinearity was extended 
for 2D discrete solitons mi, as well as for 2D systems with the spatially modulated strength of the long-range dipole- 
dipole repulsion [42]. This scheme gives rise to a variety of sturdy multidimensional states, including 2D and 3D vortex 
solitons isniEsi and more complex 3D modes, viz., soliton gyroscopes [32], vertical vortex-antivortex hybrids [33] . 
and vortex tori with intrinsic twist (“hopfions”), which carry two topological numbers [34]. These self-trapped states 
feature exceptional robustness in comparison with models of other types, as some states, such as the vortex-antivortex 
hybrids [33] and hopfions [34], do not exist in other single-component models, while multidimensional vortex solitons, 
and soliton gyroscopes, are stable only in the present setting. 

The objective of the present work is to expand the class of robust 2D and 3D modes supported by the spatially grow¬ 
ing repulsive nonlinearity, by adding other obviously important soliton species, such as the above-mentioned multipoles 
(2D dipoles and quadrupoles, and 3D octupoles), and multiplets (bound states) built of 2D vortices and antivortices, 
including vortex-antivortex pairs (VAPs) and quadruplets (VAQs). Spinning dipoles and vortex-anti vortex multi¬ 
plets are introduced too. The consideration is performed by means of numerical methods, in combination with the 
Thomas-Eermi approximation (TEA), which makes it possible to predict basic results in an analytical form. In par- 
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ticular, we demonstrate that there is a bifurcation which destabilizes 2D dipoles, replacing them by stable VAPs, 
and another bifurcation, which replaces isotropic 2D VAQs by stable anisotropic modes of the same type (AVAQs). 
Thus, the present model offers the simplest setting in which stable multipoles and vortex-antivortex multiplets can be 
constructed in the multidimensional space. In this respect, it may be compared to field-theory models, where soliton 
complexes are created in sophisticated multi-component settings [71 [35] . 


II. THE GOVERNING EQUATIONS 


The basic model is represented by the scaled D-dimensional Gross-Pitaevskii equation (GPE) which governs the 
evolution of the mean-field wave function, '0(r, t), in BEG: 


V72 / 


r(r) 


( 1 ) 


Here t is time, is the 3D or 2D Laplacian, and cr(r) represents the spatially growing local strength of the repulsive 
nonlinearity. Eollowing Ref. |30], we here adopt the steep modulation profile, cr(r) = exp(r^/2ro), where we set ro = 1 
by means of straightforward rescaling. As well as in Ref. m, milder profiles, characterized by cr(r) growing as 
(with a > D, as said above) are possible too, but the anti-Gaussian one makes it possible to present results in a 
compact form. Dynamical invariants of Eq. 0 are the norm, the angular momentum, and the Hamiltonian: 
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In the 2D reduction of the model, the 


vectorial angular momentum 0 is replaced by its single component. 
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The 2D version of Eq. 0> with time substituted by the propagation distance (z), applies as well to the spatial- 
domain evolution of the amplitude of electromagnetic waves in a bulk optical waveguide [T] with the self-defocusing 
cubic nonlinearity. In that case, the 2D reduction of Eq. 0 determines the total power of the optical beam, while 
Eq. 0 is the beam’s orbital angular momentum m- 

Stationary states with real chemical potential /i > 0 (in optics, —y is the propagation constant) are looked for in 
the usual form. 


ip (r, t) = exp {—ifit) U (r), 

where the (generally, complex) spatial wave function obeys the equation 


_ ^ 1 
^ 1 grj^2 j, Qj, ^ j,2 g^2 


U + exp(-)|UrU 


in the 2D setting, with angular coordinate cp. The 3D stationary equation is 
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where P is the usual angular part of the 3D Laplacian 
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III. MULTIPOLES: NUMERICAL RESULTS AND THE THOMAS-FERMI APPROXIMATION (TEA) 

Stationary solutions of Eq. 0 and Eq. ^ for multipole modes were obtained by means of the imaginary-time 
method mi, which is capable to generate not only the ground state but, under special conditions, higher-order modes 
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(a) (b) 


(c) 


FIG. 1. (Color online) Generic examples of numerically found dipole (a) and quadrupole (b) modes in the 2D setting, and of 
the octupole (c) in 3D. These modes are described by real wave functions U with /x = 2, see Eq. In (c), the 3D image of 
the octupole is displayed by means of red and blue isosurfaces, which correspond, respectively, to 1/ = +0.01 and U = —0.01. 



too [45], and also by means of the modified squared-operator method [46]. The respective inputs, emulating the 
multipoles sought for, were built as combinations of identical Gaussians with separated centers and alternating signs. 
Typical numerically generated examples of stationary multipoles, viz.^ 2D dipole and quadrupole, and a 3D octupole, 
are displayed in Figs. Ba), (b), and (c), respectively, for /x = 2. 

The stability of the so generated modes was examined via direct simulations of the perturbed evolution, using the 
standard pseudospectral split-step fast-Fourier-transform method. 2D simulations were performed in the domain of 
size (127r)^, covered by a mesh of 512^ points. For the 3D simulations of the octupoles, the split-step algorithm was 
used too, in the domain of size (127r)^, with the mesh of 192^ points. In particular, the course of the evolution, under 
the influence of the small numerical noise due to the discrete nature of the mesh, the octupole with /x = 2 kept its 3D 
shape virtually invariant for t = 200, which exceeds 10 respective diffraction lengths. 

Numerically found shapes of the 2D dipole {m = 1) and quadrupole {m = 2) modes, which are displayed in Figs, 
[^a) and (b), suggest that these modes may be approximated, in the simplest form, by a real ansatz^ 

Um{r, +) = V 2 B(r) cos {mcp ). (9) 

The substitution of this into Eq. (|^, and the use of the usual mean-field approximation, cos^ (mcp) (3/4) cos {mcp)^ 
yields the following radial equation: 
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The application of the TFA to Eq. (10) implies neglecting the derivative terms, which yields the simple result: 


+D(b = 


0, r2 < 


-V I , r" > ^ • 


(M-+)exp(-4), 

This approximation predicts that maxima of the amplitude are located at distance from the origin 


^(max) 

'2D 
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( 12 ) 


Amplitude maxima of the numerical solutions, found for /i = 2, are located at distances = 1.37, '^q^admpoie “ 


2.2, and r 


(max) _ 


octupole 


= 3.46 from the origin, while the TFA results (^12k and rtl7k (see below) predict, for the same case 


(m = 2), = 1.13, rl^^dlnpoie = l-8> and = 2.75. Naturally, the agreement improves for larger /u, i.e., 

larger N, as the TFA is more accurate for stronger nonlinearity (the increase of /x from 2 to 10 leads to a decrease of 
the relative difference between the numerical solutions and their TFA counterparts by a factor 2). 
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Families of the modes are characterized by the norm (total power) as a function of ja. The TFA based on Eq. 
produces the following dependence: 


»£"'(»■;-») = Tr 


- — ) e-^dp, 


(13) 


where p = r^/2, and the integral can be easily computed numerically. In the limit of /x —oo, Eq. (13) yields an 
asymptotically constant slope, 


d ^r(TFA) / , N 47r 

-]V,„ ^ oo) = 


( 14 ) 


Dependences (13) for m = 1 and m = 2 are displayed in Fig. [^a), along with their counterparts produced by full 
numerical solutions of Eq. 0- 

In the 3D geometry, the typical shape of the octupoles, displayed in Eig. [^c), suggests that the corresponding 
angular structure may be approximated by spherical harmonics with quantum numbers / = 3,m = ±2 [47]: ^32 = 
= Ccos'd sin^'d cos ( 2 (/:?), where the constant C is not essential here, 'd is the second angular coordinate, 
and Py 32 = 12132 . Thus, the ansatz for the octupole is adopted as 


U{r,'d,ip) = V 3 D {r) cos'd sin^ 'd cos {2ip ). 


(15) 


Eor the application of the mean-field approximation to the angular dependence in the cubic term of Eq. 01 
[I 32 ('^^5 ^)]^ must be projected back onto T 32 ('d, (^), which is done with the help of the standard formula, (/('d) | |^('d)) = 
fo fWdW (sin'd) d'd. The eventual results produced by the 3D version of the TEA are 
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cf. Eq. (12), and 
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cf. Eq. (13), with 
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3D 


(^^cx)) = ^(2.)3A 


(17) 


(18) 


(19) 


cf. Eq. (14). The plot produced by Eq. (18) is displayed in Eig. [^b), along with its counterpart generated by the 
full numerical solution. 

An essential finding revealed by the systematic numerical analysis of the 2D setting is that the increase of /i (and 
N) leads to destabilization of the 2D dipole mode at a critical point. 


= 1.64, = 1.38, 


( 20 ) 


by a bifurcation^ which gives rise to a new mode in the form of a YAP, which will be presented in the next section. 
Eurther, careful analysis demonstrates that both 2D quadrupoles and 3D octupoles are, strictly speaking, entirely 
unstable solutions. Nevertheless, they are robust objects, provided that their norms are relatively small, as the 
instability is very weak in that case. Eor the quadrupoles, the simulations reveal that the robustness interval is /i < 
/^max — 3, corresponding to Y < Y^ax — 1-55, see Eig. [^a). In this interval, initially perturbed quadrupoles maintain 
their shape as long as the simulations were run (with the evolution times measured in many dozens of characteristic 
dispersion times), exhibiting small-amplitude oscillations between the unperturbed shape and a coexisting stable 
isotropic VAQ mode (the latter one is presented in detail below). Eor the octupoles a similarly defined robustness 
border is located at /imax — 2, corresponding to a low value of the 3D norm, N 0.03. 
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(a) 


(b) 


FIG. 2. (Color online) Lines labeled “Dipole”, “Quadrupole”, and “Octupole” show, respectively, the 2D and 3D norms, N, 
vs. /i for these modes, as obtained from numerical computations. The lines marked by “TF” display the predictions for the 
same dependences produced by the Thomas-Fermi approximation, as per Eqs. ( |13[ ) and (|18| ). The short solid segment of the 
N{fi) curve for the dipole-mode family designates a completely stable subfamily (rf. Fig. fsTbelow). 
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(d) 


FIG. 3. Persistent rotation of a stable dipole with fi 

p = 1, = 10. 


1.5, initiated by the application of a torque, given by Eq. (21) with 


The application of the torque to the stable 2D dipole mode, i.e., the multiplication of the respective stationary 
solution by the phase factor 


T = exp (ipp tanh (x/xo)), (21) 

with real constants p and xq [32], readily sets it in persistent rotation, thus creating a robust “propeller” mode, cf. 
Ref. [T4|, which preserves the original two-peak shape. The spinning regime is illustrated by a set of snapshots in Fig. 

where the rotation period is T = 3.9. Rotation of the dipole was simulated up to t = 100, amplitude oscillations 
of the dipole in the course of rotation being limited to ^ 2% of its initial value. On the other hand, the application 
of the torque to the quadrupoles and octupoles does not generate persistently rotating states, in accordance with the 
above conclusion that these species do not represent fully stable modes. 


IV. VORTEX-ANTIVORTEX PAIRS (VAPS) 

In addition to multipoles, bound pairs of vortices with opposite signs of angular momenta were investigated. A 
typical example of the VAP and the respective bifurcation diagram are displayed in Figs. |^and[^ respectively. As 
seen in panel (a) of Fig. the amplitude profile of the VAP is quite similar to that of the dipole mode, but panel (b) of 
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|U| Phase{U} 



(a) (b) 


FIG. 4. An example of the amplitude and phase profiles of a 2D stable vortex-anti vortex pair (YAP), obtained in the numerical 
form for /x = 5 (the corresponding norm is N = 14.78). The phase distribution is displayed in degrees. 



FIG. 5. N(fi) curves for the dipole-mode and vortex-antivortex-pair (YAP) families. Stable and 
indicated by continuous and dotted lines, respectively. The red dot marks the bifurcation point (20), 
branches off from the dipole-mode one. 


unstable subfamilies are 
at which the YAP family 


the same figure demonstrates the phase structure of the stationary YAP, which its real dipole-mode counterpart does 
not have. It is exactly the phase structure which makes it possible to identify the new mode as the vortex-antivortex 
bound state. 


In accordance with what is said above, the stability of the modes under the consideration was identified by means 
of systematic simulations of their perturbed evolution. As shown in Fig. the YAP is stable after the bifurcation 
point, viz.^ at 1.64 < fi < 4.06, 1.38 < N < 10.35, then it is destabilized by oscillatory (in t) perturbations in a 
finite region adjacent to the stability area, 4.06 < /i < 8.06, 10.35 < N < 30.76, which is followed by the eventual 
restabilization of the YAP at /i > 8.06, N > 30.76. Stable YAPs can be easily set in spinning motion, similar to what 
is shown above for dipoles. 


The emergence of the YAPs determines the instability of the dipole modes. In the interval adjacent to the critical 
point Eq. (20), 1.64 < /x < 3, the dipole features a weak instability, periodically oscillating between its unperturbed 
shape and the coexisting stable YAP, whose amplitude profile is close to the dipole’s one (not shown in detail here). 
A similar dynamical regime, featuring remittent shape revivals of a weakly unstable dipole state, was observed in a 
model with a nonlocal nonlinearity m- At /X > 3, the dipole’s instability gets stronger, leading to its spontaneous 
transformation into a fundamental isotropic state (not shown in detail either). 











|U| Phase{U} 



(a) (b) 


FIG. 6. The amplitude (a) and phase (b) profiles of a stable isotropic vortex-antovortex quadruplet (VAQ) at /x = 5 and 
N = 8.12. 





FIG. 7. Juxtaposed N(fi) curves for the (fully unstable) quadrupoles, as produced by the numerical results and by the Thomas- 
Fermi approximation (“TF”) [these two curves are the same as in Fig. [^a)], and for the (partly stable) isotropic and (fully 
unstable) anisotropic vortex-antivortex quadruplets (VAQs). As usual, stable and unstable families of solutions are indicated 
by solid and dotted lines, respectively. The dashed-dotted and dashed line s refer, severally, to the Thomas-Fermi approximation 
for the quadrupoles and the anisotropic VAQs. The bifurcation point (22), at which the anisotropic-VAQ branch splits off from 
its isotropic counterpart, is marked by the red dot. 


V. VORTEX-ANTIVORTEX QUADRUPLETS (VAQS) 


Similar to the 2D dipoles, which coexist with VAPs, real stationary 2D quadrupoles coexist with a branch of stable 
complex solutions for vortex-anti vortex bound states in the form of the {isotropic) VAQ, see a typical example of the 
latter in Fig. However, Fig. [^demonstrates the difference from the situation for the coexistence of the dipoles and 
VAPs: the isotropic-VAQ branch emerges at /i = 0, rather than branching off from the quadrupole one at a finite 
value of /X, cf. Fig. Recall that, unlike the dipoles, quadrupole solutions are, strictly speaking, always unstable, 
hence, indeed, a new stable branch cannot bifurcate from them. As shown in Fig. [^ , the isotropic-VAQ family is 
stable at 0 < /X < 5.74 (which corresponds to 0 < V < 11.06). At /x > 5.74 {N > 11.06), these modes become unstable 
against oscillatory perturbations, eventually evolving into the isotropic ground state (not shown here in detail). 

Further, a new branch, of anisotropic VAQs {AVAQs)^ bifurcates from the family of their isotropic counterparts at 


^(aniso) ^ j^ianiso) ^ q_33_ 


( 22 ) 


A typical example of an AVAQ mode is displayed in Fig. The analysis demonstrates that this solution is always 
unstable [this fact explains why the bifurcation occurring at point (22) does not destabilize the branch of the isotropic 
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(a) (b) 


FIG. 8. The same as in Fig. but for an unstable anisotropic VAQ, with // = 5, iV = 7.40. 
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FIG. 9. The transient rotation of the unstable AVAQ, with fi — 12, N — 38.08, prior to its spontaneous transformation into a 
stable vortex-antivortex pair. 


VAQs]. Specifically, from the bifurcation point (22) up to /i = 12.54 (A^ = 40.77), a perturbed AVAQ solution starts 
spinning motion, with an angular velocity determined by the initial perturbation, see an example in Fig. (in this 
case, the conservation of the angular momentum is provided by the recoil effect of small-amplitude waves emitted by 
the perturbed AVAQ, which are not visible in Fig. [^. Eventually, the unstable spinning AVAQ transforms not into 
the isotropic ground state, which is typical for other species of unstable modes which were mentioned above, but, 
instead, into a stably rotating VAP. In several cases that were examined, this transition occurred before the AVAQ 
would complete half a cycle of its rotation. Finally, at /i > 12.54 (V > 40.77), the AVAQ immediately transforms 
into a rotating VAP, without going through the intermediate spinning stage (not shown here in detail). 


VI. CONCLUSION 

The recently introduced class of optical and matter-wave models, with the strength of the self-repulsive cubic 
nonlinearity growing from the center to periphery, gives rise to a large variety of completely stable complex 2D and 
3D self-trapped states (solitons), which do not exist or are strongly unstable in other physical settings. Here we have 
demonstrated that, together with previously found states, the same models support several other species of 2D and 
3D solitons which are not available (in a stable form) in other systems either. These are 2D dipoles and quadrupoles 
and 3D octupoles, as well as VAPs (vortex-antivortex pairs) and isotropic or anisotropic VAQs (vortex-antivortex 
quadruplets). The modes found here are stable or weakly unstable, surviving for a long time (over a long propagation 
distance), which makes them physically relevant objects. In addition to quiescent states, persistently spinning dipoles 
have been found too. The results, which were obtained by means of numerical methods and analytically, using the 
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TFA (Thomas-Fermi approximation), demonstrate the potential of the settings with the spatially growing strength 
of the self-repulsion for the creation of a great variety of complex stable multidimensional modes. 

Natural extensions of the present analysis may be additional analysis of the 3D setting (in particular, for 3D dipoles 
and quadrupoles), and its development for two-component systems, in which the multitude of nontrivial self-trapped 
states may be further expanded, as suggested, e.g., by results for the two-component system in ID [48] , 
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